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Abstract. This paper studies the problem of exact localization of multiple objects with 
noisy data. The crux of the proposed approach consists of random illumination. Two 
recovery methods are analyzed: the Lasso and the One-Step Thresholding (OST). 

For independent random probes, it is shown that both recovery methods can localize 
exactly s — 0{m), up to a logarithmic factor, objects where m is the number of data. 
Moreover, when the number of random probes is large the Lasso with random illumination 
has a performance guarantee for superrcsolution, beating the Rayleigh resolution limit. 
Numerical evidence confirms the predictions and indicates that the performance of the 
Lasso is superior to that of the OST for the proposed set-up with random illumination. 



Two-point resolution is a standard criterion for evaluation of imaging systems, i.e. the 
ability of the imaging system to distinguish two closely located point objects. The smallest 
resolvable distance £ between two objects, called the (two-point) resolution length, is then 
defined as a metric of the resolving power of the imaging system. Let A be the aperture 
of the imaging system, z the distance to the objects and A the wavelength. The classical 
Rayleigh resolution criterion then states 



where there is some arbitrariness in the constant depending on the precise definition of 
minimum resolvable length £. 

For noisy data, such a criterion is more difficult to apply as determination of I becomes a 
statistical problem. One option would be to formulate the two-point resolution problem as 
a statistical- hypothesis-testing problem (one versus two objects), see [23J, [33] and references 
therein. However, it is cumbersome to generalize this approach to multiple point objects. 

In this paper we first study the resolution issue from the perspective of exact, simultane- 
ous localization of multiple point objects. We evaluate an imaging method by saying that 
it can exactly localize s (sparsity) randomly distributed point objects mutually separated 
by a minimum distance I with high probability. In addition to reconsidering the issue of 
resolution, we seek an approach that can recover a high number s = 0(m) objects where m 
is the number of data, with resolution I far below what is dictated by the Rayleigh resolution 
limit Q (see Remark [3]). This latter effect is called superresolution. 

Consider the noisy data model: 
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where X G C N is the object to be recovered, Y G C m is the data vector and E G C N 
represents noise. We shall assume that <fr has unit-norm columns. This can always be 
realized by redefining the object vector X. 

Sparse object reconstruction for this model can be broken into two steps: localization 
(i.e. support recovery) and strength estimation. For underdetermined sytems, the former, 
being combinatorial in nature, is by far more difficult than the latter which is a straight- 
forward inversion if the former is exact. The former step is called model- selection in linear 
regression and machine learning theory [3j El HUJ [121 EH [28j [351 E] from which one of the 
reconstruction methods studied in the present paper originates. 

Exact localization with noisy data is challenging. Many reconstruction methods guarantee 
stability (i.e. the reconstruction error bounded by a constant multiple of the noise level) but 
not necessarily exact localization. Orthogonal Matching Pursuit (OMP) is a simple greedy 
algorithm with proven guarantee of exact localization for sufficiently small noise and worst- 
case coherence. 

A basic quantity for stability analysis in compressed sensing is the notion of coherence. 
Let the worst-case coherence be defined as 

I $*$,-! 

(3) M*) = max J 



i+3 || *j || 2 || || 2 
A standard result is the following result [18] . 



Proposition 1. Consider the signal model |5|j. Suppose the sparsity s of the real-valued 
object vector X G IR satisfies 

lie 

s<-(l + -) — , X min = mm\Xi\ 

A [i /iA min ies 

Denote by X s the output of OMP which stops as soon as the residual error (in £ 2 -norm) is 
no greater than e. Then 

(i) X s has the correct support, i.e. 

supp(A £ ) = supp(A) 

(ii) X s approximates the true object vector 

\\X £ -X\\l< 



12 ^ 1 -/i(s- 1) 
The general lower bound [T71 HO] 



N — m 



m(N- 1) 

for the mutual coherence of any m x N matrix $ implies that the sparsity s allowed by 
Proposition [l] is 0(\fm) for N 3> m. 

A main purpose of the paper is to explore the utility of two other methods from compressed 
sensing theory, the One-Step Thresholding (OST) [3j and the Lasso [351 CS], that have the 
potential for exact localization of much higher number 0(m) of objects. 
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The One-Step Thresholding (OST), proposed in [3], involves just one matrix multiplication 
plus thresholding: Compute Z = and determine the set of points 

S = {i e {l, ...,#} : \Zi\>n} 

for some threshold r*. In other words, the OST is the linear processor of Matched Field 
Processing (MFP) plus a thresholding step [2]. On the other hand the linear processor of 
MFP is the same as the first iterate of OMP. Consequently OST has even lower complexity 
than OMP which is its main appeal. 

For the OST's performance guarantee, we need the notion of average coherence defined as 

i 



= max 

1 ; N-l f 



E*^ 

fry 



in addition to the worst-case coherence. 

The following is the performance guarantee for OST [3]. 

Proposition 2. Consider the signal model Assume that X £ M. N is drawn from the 
generic s-sparse ensemble of real-valued objects. Assume E to be distributed as CN(0,a 2 I) ; 
the complex Gaussian random vectors with the covariance matrix a 2 l. 
Suppose 

(4) < ' 



fm ~ V10 log N 
for some C\ > ( which may depend on log N ) and 

(5) *(•) < 12 " ( * ) 



r m 

Assume \\X\\2 = 1. Define the threshold 

(6) r* = 4 v / loi max | a, 

Suppose the number of objects obeying 

m 

and i/iai 

(8) X min = min \Xi\ > 2r*. 

T/ien i/ie OST threshold satisfies F (s ^ s\ <9/N. 

In other words, for sufficiently small worst-case coherence Q and average coherence ([5]) 
and noise dSl) , OST can exactly localize C(m) objects, up to a logarithmic factor with high 
probability. Once the support is exactly recovered, an estimate X can be obtained by 
pseudo-inversion on the object support S. 

The other method studied in this paper is the Lasso [35]. The Lasso estimate X is defined 
as the solution to 

(9) min h\Y - $Z\\ 2 2 +<ya\\Z\\ u 7>0 

Z Ji 
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where 7 is a regularization parameter. 

The following sufficient condition for exact localization by the Lasso is given by 

Proposition 3. Consider the signal model Assume that X G M. N is drawn from the 
generic s-sparse ensemble of real-valued objects. Assume E to be distributed as CN(0,a 2 I). 
Suppose that <fr obeys the coherence property 

< 10 > 5 if* 

with some positive constant clq. Suppose 
(11) s< <*" 



1*11! log jv 

for some positive constant Cq. Let S be the support of X and suppose 

(12) X min > 8^2 log X. 
Then the Lasso estimate X with 7 = 2 a/2 log N obeys 

(13) supp(X) = supp(X) 

(14) sign(X) = sign(X) 

with probability at least 1 - 2X- 1 ((2tt \ogNy 1 ' 2 + sN' 1 ) - £>(X- 21 °s 2 )). 

Some comparison between Proposition [3] and [2] is in order. Both deal with randomly 



distributed objects. Both J4J) and (10) are sufficiently weak assumptions for most imaging 
problems. Also (12) and (pF are similar when /1 = 0(a). The lower bounds for the success 
probabilities are comparable up to a logarithmic factor. The main technical assumption of 
Proposition [3] is ( fTT| ) while for Proposition [2] it is ([5]). When the operator norm ||*||2 obeys 
the bound ||*||2 = 0(N/m) condition ( fTT| ) is comparable to ([7]). 

A drawback to Proposition's that the thresholding rule (|6j) requires the precise knowledge 
of \i which can only be calculated numerically. As we shall see, the Lasso-based method also 
has a better numerical performance than does the OST (cf. Figure [5] and [7]). 

To realize the potential of the two above results in imaging, we shall consider the idea of 
random illumination for point scatterers. We shall show that a suitable condition of random 
illumination enables us to (i) obtain a guaranteed exact localization of s = 0(m), up to a 
logarithmic factor, objects and to (ii) harness the superresolution capability (i.e. breaking 
the Rayleigh resolution limit ([!])). 

Previously we have studied the problems of imaging point scatterers [21] [23] using coher- 
ence and operator-norm bounds. We shall demonstrate that the imaging performance can be 
significantly improved by random illumination. In particular, suitable random illumination 
leads to superresolution. 

However, both Propositions [2] and [3] share the following common drawbacks: (i) they are 
restricted to random objects; (ii) they do not address the reconstruction error when the error 
level is above threshold and exact localization is unattainable; (iii) they are limited to the 
i.i.d. Gaussian noise model. Issue (i) is pertinent particularly to imaging extended objects 
whose supports are clearly not random. Issue (ii) is related to robustness with respect to 
a wider range of error. Issue (iii) arises in optics where the Poisson or shot noise model is 
more appropriate. 
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The standard compressed sensing method that are without any of the above limitations 
is the Basis Pursuit Denoising (BPDN) 



(15) min||Z||i, s.t. \\Y - *Z|| 2 < £ 

z 



[To] . BPDN, of course, is equivalent to the Lasso (55) for an appropriately chosen 7. 

The performance guarantee for BPDN is typically given in terms of the restricted isometry 
property (RIP) due to Candes and Tao [13]. Precisely, let the sparsity s of a vector Z G C N 
be the number of nonzero components of Z and define the restricted isometry constant (RIC) 
5 S G [0, 1] to be the smallest nonnegative number such that the inequality 

(16) (l-6 s )\\Z\\t<\\*Z\\l<(l + 6 s )\\Z\\l 

holds for all Z G of sparsity at most s. BPDN has the following performance guarantee 

[Uj. 

Proposition 4. Suppose the RIC satisfies the bound 

(17) 5 2s <V2-l. 

Then the BPDN minimizer X is unique and satisfies the error bound 

\\X - X\\ 2 < C lS - 1/2 ||X - + C 2 e 

where X^ s > is the best s-sparse approximation of X and Ci,C 2 are absolute constants de- 
pending on 5 2s only. 

In Proposition [4| BPDN does not guarantee the exact recovery of the discrete support X, 
which is less important for extended objects, but also does not have any of the limitations 
mentioned above for Propositions [2] and [3j 

The plan for the rest of the paper is as follows. In Section [2j we review the forward 
scattering problem and the paraxial approximation. We describe the set-up of random 
illumination in the paraxial regime. In Section [3] we state and prove the main results. 
In Section [4] we analyze the performance of BPDN with random illumination for extended 
objects in Section[4]and discuss the issue of resolution in imaging extended objects. In Section 
[5] we give the worst-case coherence bounds. In Section [6j we give the average coherence 
bound. In Section [7] we give an operator norm bound of 0(N/m) required to guarantee a 
nearly optimal performance for the Lasso. In Section [8] we present numerical simulations to 
verify the predictions and show the superiority of the Lasso over the OST for the set-up of 
random illumination. We also present numerical results for extended objects. We conclude 
in Section [9l 

2. Point scatterers and paraxial approximations 
Let £ be a finite square lattice of spacing £ in the object plane {z = 0} C M 3 : 

(18) C = {r ; : I = 1, N} = :i,j = 1, Tiv} , I = (i - 1)VN + j 

and suppose that s point scatterers are located at grid points of C. The total number of grid 
points N is a perfect square. 
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Let Tj 6 C, I — 1, N be the reflectivity of the scatterers. The scattered field u s obeys 

N 

(19) u s (r) = ^TjGfaTjWiTj) +u s (v j )) 
for any r {r^ : 7^ 0} where -u 1 is the incident field and 

iw|r-r'| 

(20) G(r, r') = — r, Vr ^ r' G M 3 

47r|r — r'| 

is the Green function of the operator — (A + cu 2 ). 

In the Born scattering approximation, u s on the right hand side of (19) is neglected, 
resulting in 

N 

(21) M s (r) = ^r,G(r,r,K(r,) 

Let aj,j = 1, ...,n be the locations of the sensors in the sensor plane {z = zq} C IR 3 
and write a,,- = (£,j,r]j, zq) where £j and rjj are chosen independently and uniformly from the 
discrete subset of [0, A] 

(22) V={-^:q=l,...,VN 



'N 

where A is the aperture of the sensor array. 

In the Fresnel approximation under the condition 



(23) ^ + T> 4 « 1 

\z$ 

the Green function G can be approximated by 

(24) G par (r,a) = e H-SI 2 /(2^) e H^l7(2*o) ; r = ( ), a = (^^zq), 

47T2 

called the paraxial Green function. 

In the subsequent analysis we shall assume both the Born and paraxial approximations in 
the scattering model. 

A main ingredient of the proposed approach is random illumination which has recently 
been used extensively for wavefront reconstruction and imaging [U [13 ED] • Here we consider 
random phase modulation (RPM) which is a random perturbation of the phase of a wavefront 
while maintaining the amplitude of the near field beam almost constant. The advantage of 
phase modulation, compared to amplitude modulation, is the lossless energy transmission 
of an incident wavefront through the modulator. In optics RPM can be created by random 
phase plates, digital holograms or liquid crystal panels [H |3l] . 

3. Main results for point objects 

We assume that as a result of p independent realizations of random phase modulators the 
incident field at the grid points can be represented as e k i,k = l,...,p,j = 1,...,N where 
9kj are i.i.d uniform random variables in [0, 2tt] (i.e. circularly symmetric). The information 
about 9kj is incorporated in the sensing matrix. 
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Figure 1. The imaging geometry for point objects 



Let the scattered field u\ is measured and collected by n sensors located at aj, I = 1, n. 
Let X = (r 3 -)f e C N be the object vector and Y = (ft) = (u s k (a.j)) G C np , i = (k-l)n+j,j = 
1, ...,n, the data vector. 

After proper normalization, the data vector Y can be written as ^ with the sensing 
matrix $ being the column-normalized version of [G paT (a.i,rj)u k (rj)], i.e. 

(25) <fnj = -J— e ^\xj-(i\ 2 /(^o) e iu,\ yj -rn\ 2 /(2z ) e ie k ^ i = ( k _^ n + L 



np 

Here m = np is the number of data. 

Our first result is a performance guarantee for the OST with random illumination in the 
diffraction-limited case satisfying the Rayleigh resolution criterion. 



Theorem 1. Let 

(26) 

Suppose 

(27) 

and 

(28) 

Then with probability at least 



N 2 <°-e K2 l\ 5,K>0. 



np > 40K 4 log N 



M_ 



1. 



(29) 



1-25-4* 



OST with the threshold 



7T 

can 



4 

VP 



8Ne 



-12t 



2 / jy-i 

np 



V* > 



localize exactly s objects satisfying 



Remark 1. The constants 5 and K in (26) are controlling parameters. 5 can be adjusted 



to control the lower bound (29) for success probability and then K can be adjusted to control 



the number of grid points in the computation domain and the number of data. 



For example, suppose 5 = 1% is acceptable. Then (26) with K = 10 implies a computation 
domain of about 0.1e 25 /\/2 grid points. 
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Proof. The proof of Theorem [T] relies on Proposition [2] and the following three lemmas. 
Lemma 1. Under (26), the worst case coherence satisfies 

'2 



P 



aKy/2 2K 



> 1 -25 



(30) _ lf 

y'p ^/np 

where a is given by (35). 

In particular, if (28) holds then a = and (Sty becomes 

(31) P {//(*) < 2K 2 / y /rip} > 1 - 25. 

The proof of Lemma[l]is given in Section|5j The utility of estimate (30 ) lies in the situation 
where both the aperture and the sensor number are limited but the number of probe waves 
is exceedingly large (see Remark |3|. For the proof of Theorem [T] we need the estimate (31 ). 

Lemma 2. 

Lemma [2] is an easy consequence of the Berry-Esseen theorem and its proof is given in 
Section 15.21 



(32) 



P 



Under the assumption (28), 
2t x t 2 



M*)> 



> l — 2ti 



2t, 



n 



Lemma 3. Let (28) hold true. Then for any c > 

c 



(33) 



P<^ v(9) < 



np 



> 1 - 8iVe" 



' N-l 
np 



The proof of Lemma [3] is given in Section [6} 

First of all, by the upper bound (30) for the worst case coherence and setting c\ = 2K 2 
in Q then the first inequality of (|4j) holds with probability at least 1 — 25. The second 
inequality of ^ follows from (26) and (27) and holds with probability at least 1 — 25. 

Second, the lower bound (32) for the worst case coherence, with t\ — t% — t, and the upper 
bound (33), with c = 24t 2 , for the average coherence imply that (|5| holds with probability 
at least 

^ 4 4 



I -At 



8Ne 



-12^ 



I N-l 
np 



n 



This completes the proof of Theorem [T] □ 
Our second result is a performance guarantee for the Lasso with random illumination. 



Theorem 2. Let (26) hold and suppose 

aKy/2 



(34) 

where 
(35) 



2K 2 
+ < 



a 



np log iV 



max IE ( e ^M^'-^)/M E UvMv, 



-Vj)/zo^ 



Assume that the s objects are real-valued and satisfy (12) and 

c np 



(36) 



s < 



21ogiV' 

8 



Then the Lasso estimate X with 7 
at least 



2y/T\ogN has the same support as X with probability 



(37) 



1 —25 — pn(n — 1 



7T 



np — 1 



2n 2 p(p 



2 V N 

-2^ 1 ((2vrlogiV)^ 1/2 + sN^ 1 ) - 0(iV~ 2 



\\ e (np-l)2 
2 ))- 



Remark 2. While it requires that N np for the bound (2£ty to approach unity, it demands 
a much stronger assumption N 3> max{pn 5 , p 2 n 2 } for the bound (31) to behave the same 
way. Numerical evidences indicate the latter to be a pessimistic estimate. 

For the special case of single sensor n = 1, the probability lower bound (31) is substantially 
improved and requires N 3> p 2 to approach unity. On the other hand, for 12$ to approach 
unity, it is necessary that n — > 00 (hence n = 1 is not an option). 

Remark 3. The superresolution effect can occur when the number p of random probes is 
large. Consider, for example, the case of n = 1 and hence the aperture A is essentially zero. 
Since a < 1, the condition 



Ky/2 + 2K 2 



and 



VP 



s < 



< 



«0 



logiV 



cop 



implies that the Lasso with 7 
probability at least that given by (31) 



21ogiV 

2y/2 logiV recovers exactly the support of s objects with 



This superresolution effect should be compared to that with deterministic near-field illumi- 
nation [20] . 

Proof. The proof of Theorem [2] uses Proposition [3j Lemma[T]and the following operator-norm 
bound. 



Lemma 4. We have 



(38) 



P 



1*112 < 



2N 
np 



> 1 — pnin — 1) 



TlyJ np — 1 



2ViV 



2n 2 p(p — l)e 



N 

(np-1) 2 



On one hand, Lemma [T] and (34) imply that (10) holds with probability at least 1 — 25. 
On the other hand, Lemma 4 and ((36|) imply that (11) holds with probability at least 



given by the right hand side of ( 38 ) . 



2. 



Combining the two and using Proposition [3] we obtain the desired statement of Theorem 

□ 



To further demonstrate the advantage of random illumination, let us consider the imaging 
set-up of multistatic responses (MR) which consists of an array of n fixed transceivers which 
are both sources and sensors (i.e. transceivers). One by one, each transceiver of the array 
emits an impulse and the entire array of transceivers records the echo. Each transmitter- 
receiver pair gives rise to a datum and there are altogether n 2 data forming a data matrix 
called the multistatic response matrix. By the reciprocity of the wave equation, the MR 
matrix is symmetric and hence has at most n(n + l)/2 degrees of freedom. 
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Recalling the coherence and operator norm bounds established in [21] and using Proposi- 
tion [3] as in the proof of Theorem [2] (below), we have the following result [21] analogous to 
Theorem |2] 

Proposition 5. Let the locations of the n transceivers be i.i.d. uniform random variables 
in [0,A] 2 . Let Q and (g) hold true. 
Suppose 

K 2 log N 



n > 



a 



and that the s real-valued objects satisfy ( |1£| ) and 

c n(n + 1) 



(39) 



s < 



4\ogN 



Then the Lasso estimate X with 7 = 2^2 log has the same support as X with probability 
at least 



(40) 1-2V25 



pn 5 / 2 (n + I) 5 / 2 

^25/2^1/2 



2N-\(2nlogN)- 1 / 2 + sN' 1 ) - £>(A^ 21og2 )). 



Remark 4. The main drawback of the lower bound ^es in the third term which requires 
N ^> n 10 to diminish. 

More generally, one can consider the case of p transmitters and n receivers, all randomly 
and independently distributed in [0,A] 2 . Then an extension of the bound (40), which is 
omitted here, requires N 3> n 5 p 5 (cf. |21j ). 



From dimension count, a fair comparison with Theorem |2| would be to set p = (n + 
l)/2 and match their degrees of freedom, i.e. n{n + l)/2. However, Proposition [5] does 
not guarantee superresolution when (28) is violated preventing the worst case coherence 
from being sufficiently small due to the deterministic nature of the illumination. Also, the 
probability lower bound (40) has a less favorable scaling behavior (N ^> n 10 ) than ([37]) for 
p = (n+ l)/2 (N ^> n 6 , cf. Remark [2]). Indeed, the numerical simulations show the recovery 
with random illumination has a higher success rate than the MR recovery (Figures [3] and 111) . 



4. Sparse extended objects 

We extend the above results to the case of sparse extended objects here (Figure [2]). 
We pixelate the sparse extended object with N pixels — 1, N of size £ to create a 
piecewise constant approximation of the object. The centers of the pixels are identified as 



£ given in (18). Let O(r) be a the original object and Oi its ^-discretization, i.e. 

N 

0, = ^I D .0( ri ) 

3=1 

where In . is the indicator function of the pixel We reconstruct the discrete approximation 
Oi by determining the object function restricted to £, denoted still by X = (0(rj)), by 
compressed sensing techniques. 
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Figure 2. The imaging geometry for extended objects 

Under the random illumination u\, pixel Qj now produces a signal at the sensor a; of the 
form 

0{v 3 ) [ G par (r,a0e-^ 2 /( 2zo )e-^ 2 /( 2z °)4(x,y)^ 

where the quadratic phase factors are due to the presence of a parabolic lens immediately 
after the object plane (Figure [2]). This lens is introduced here to simplify our analysis. In 
practice, the lens is not needed and should have a negligible effect on performance. 

As for the case of point objects we assume that as a result of the RPM u\ takes a constant 
value e *•»' in pixel Qj and that 9^ are i.i.d. random variables in [0, 2tt] as a result of random 
phase modulation. 

The total signal produced by Og and detected at sensor a; is 

Vofa)^' / G P Ur,zi)e- tuJx2/(2zo) e- lujy2/{2zo) dxdy 

= Q(rj)e l8kj e lu} ^' ^ 2zq ^ e %UJ ^ ^ 2zo ^ e~ x 3 e ~ tuJVl Vj I z ° / e~ tuJ ^ lX ^ zo e~ %U)my ^ zo dxdy 

□ 



plus an error term which includes the discretization error and external noise where □ 
denotes the square of size £ centered at the origin. Since 

<«) L — — ^ - % - (f ) ^ - g) - »w 

independent of the pixel index, we can normalize the data by dividing the signal at sensor I 
by this number as long as 

Dividing the data further by the phase factors e tU) ^ /( 2z o) e tuJ Vi /(2«o) anc [ y/r^ we write the 
signal model as ([2]) with the sensing matrix element 

1 



(43) <pij = e *ki er ™ti*i/«> e -™nvi/*> i i = (k -l)n + l 



np 



n 



The difference between the signals produced by O and its discretization Og is the dis- 
cretization error -Edisc- How small must i be in order for the ^ 2 - n o rm of the discretization 
error -Edisc be less than, say, e after rewriting the signal model as ([2])? This can be estimated 
as follows. 

First, by the inequality H-Edisclh < || iodise ||oo\/^P it suffices to show ||£?disc||oo < e l ' 
Since 

N 

t&r) (r) 

3=1 

is the illumination field, the uncontaminated signal detected by sensor a; in the absence of 
external noise in the signal model ^ is 

1 N r 

(44) (JFO)i = — - V / O(x,y)e- luj ^ /zo e- lujrny ^ zo dxdy, 
for i = (k — l)n + I. On the other hand we have 

i N r 

(FO e )i = -^^^O^vV]) J e-^m/z° e -™r,m/zo dxdy 
j=i j 

for i — (k — l)n + I. By definition 

£disc =TO- TO, e C pn 

and hence 

(45) n^sciioo < l|Q " Q ; ll t; 

mm/ \g{a.i)\ 

where || • H^i denotes 

II/IU* = J \f{x,y)\dxdy, 
i.e. the norm of the function space L 1 . Therefore we have the following statement. 
Lemma 5. If 

(46) \\0-O e \\ L i < ^ minima/) | 

■Jnp i 



then 

-Edise 2 < £• 



Remark 5. The presence of the factor (np) in (4-6) is due to the transition from L 1 
function space norm to the discrete t^-norm. 



Since the sensing matrix (25) for the point objects can be written as 

Di$D 2 



where $ is as ( 43 ) and 

D x = diag(e^ 2/(2 ^ o) e^ 2/(22o) ) 
D 2 = diag(e^ 2/(22o) e^ 2/(2z(,) ^ 
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are diagonal, unitary matrices. All the preceding results, including Theorems [T] and |2j can 
be proved for the sensing matrix (43) by minor modification of the previous arguments. 

However, the object vector X = (0(i\,)) of an extended object generally does not fall 
into the category of random point objects assumed in either Proposition [2] or [3] since by 
definition the discrete approximation of an extended object must cluster in aggregates and 
its amplitude typically changes continuously. So we take an alternative approach below by 



resorting to the minimization principle (15) of BPDN. 



The RIC for a structured sensing matrix such as (43) is difficult to estimate directly except 
for the case of single shot (p = 1) and the case of one sensor (n = 1). For the one-sensor 
case, (43) with (£2,7^) = (0,0) is the complex-value version of the random i.i.d. Bernoulli 
matrix: 



(47) 



whose RIC can be easily estimated by the same argument given in [1]. The single sensor 
imaging set-up resembles that of Rice's single-pixel camera [19] which employs a discrete 
random screen instead of a random phase modulator. 



For the single-shot case, the sensing matrix (43) is equivalent to the random partial Fourier 



matrix, modulo an unitary diagonal matrix, and the standard RIP estimate [29J requires the 
Rayleigh criterion (28) to be met which guarantees (42) with probability one. However, 
there exists a small probability of a; = (£/, rji) falling near the boundary of the aperture and 
hence a small value of |g(a;)|. Normalizing the data by (^(a^l then carries a small risk of 
magnifying the errors. 

For the general set-up with multiple shots and sensors, we use the mutual coherence to 
bound the RIC trivially as follows. 

Proposition 6. For any seNwe have 

S s <^)(s-1). 

Combining Lemma [T], Propositions [4] and [6] we obtain the following result. 



the sensing matrix (43) and sparsity up to 

1 



Theorem 3. Under (26]), the RIC bound (11) holds true with probability at least 1 — 25 for 



(48) 



where 



s < 




aKy/2 2K 2 
+ 



-1 



y/P 



np 



max |E ( e *io>(*j'-*i)l*o) E ( e «»«(Vj'-w)/*> 
itt' 1 v ' v 



c.f. (S§. 

Furthermore, suppose the total error in the data is E = -Edisc + E ext where -Edisc o,nd E ext 
are, respectively, the discretization error and the external noise. Then the reconstruction X 
by BPDN satisfies the error bound 



(49) 



\X - X\\ 2 < C lS -y 2 \\X - X^\\, + C 2 (\\E disc \\ 2 + \\E cxt \\ 2 ) 



for all s satisfying (48). 
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Remark 6. Since BPDN does not guarantee exact localization, an appropriate metric for 
resolution can be formulated in terms of the smallest pixel size £ m i n and largest sparsity s 



such that (4-9) holds true with both the discretization error E^ sc and s~ l l 2 \\X — X^\\i being 
reasonably small. 

The right definition of "small errors", however, is problem specific. The discrete norms 
(t\- or £ 2 - norm) tend to go up simply because the effective sparsity increases. Hence the 
right metric of reconstruction error should be properly normalized by the size of the object. 



For example, consider the special case when X is s-sparse. Then we can rewrite (4-9) as 
(50) s- l / 2 \\X - X\\ 2 < C 2 ||£ext||2S~ 1/2 + C 2 \\E disc \\ 2 s- 1/2 

whose left hand side is a measure of the reconstruction error per pixel of size I. 

Below the diffraction limit A£/(Xzq) < 1 (a ^ 0), one can reduce the discretization error 
by reducing the pixel size according to Lemma^ On the other hand, the sparsity s increases 



in proportion to £ 2 for a two-dimensional extended object. To satisfy (48) the smallest 
admissible pixel bounded from below roughly by 

(51) £ min Z a 1 ^- 1 / 4 

meaning that the minimum super-resolved scale decreases at least as fast as the negative 
quarter power of the number of random illuminations. 
For the diffraction-limited case a = 0, we have instead 

4nn * n" 1 V 1/4 



which is more favorable than (51) for n 1. However, the discretization error bound 
(Lemma^) is less useful in this case. 

5. Worst-case coherence bound 
5.1. Proof of Lemma [if upper bound. 

Proof. Summing over a/, I = 1, n we obtain 

p n 1 p n 

k=l 1=1 U P k=l 1=1 

We shall estimate the two summations separately. 

First consider the summation over random illuminations k = 1, p. Define the random 
variables Ai, Bi, I — 1, n, as 

(53) Ai = cos [9 kj - 9 kj/ ] 

(54) Bi = sin [0 fci - 0jy/] 

and let 



(55) S p = ^2(A i + iB l ). 



i=i 



To estimate S p , we recall the Hoeffding inequality 
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Proposition 7. Let A\ + iB 1 , ...,A p + iB p be independent random variables. Assume that 
A h B[ e [ai,bi\,l = l,-..,p almost surely. Then we have 

p 2 t 2 



(56) W[\S P - E5 P | > pt] < 4 exp 

for all positive values of t. 



n - an 



We apply the Hoeffding inequality to S p with ai = —1,6/ = 1,/ = l,-..,p and 



t = Kj-, K>0 
p 



to obtain 



(57) 


P 


p- 1 


p 


>-«4 








k=l 







\S p {9kj'—0kj, 



Note the dependence of S p on 6kj — Qkj' and the symmetry: \S p (9kj> — 0kj) 
As a consequence, there may be N(N — l)/2 different values of S p . By union bound with 
(57), we obtain 



(5? 



P 



p max 



fe=i 



> K 



< 2N(N - l)e~ K2/2 < 5 



by (26). 



Next consider the summation, denoted by T n , over the sensor locations / = 1, n in (52 ): 

n 
1=1 



By the same argument we obtain 
P 

and hence 
(59) P 



maxn 1 \T n — ETJ > K 



max — \T n \ > a + K\ — 
j'¥=j n V n 



< 2N{N - l)e- R2/2 



< 8, a = max — I ETJ 
n 



by (26). 



By the mutual independence of £/ and rji we have 



a = max — IETJ 
oft' n 



1 

max — 



i=l 

ix |E (e^^^j' - ^^ 20 ) E (e^^^' - *^ 20 ) | 
since £i,r)i,l = 1, n are independently identically distributed. 



max 



Combining (59) and (58) and noting the independence of these two events, we obtain 

aKV2 2K 2 
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rap 



with probability at least 1 — 26. 

Simple calculation with the uniform distribution on the set D given in ( 22 ) yields 



(60) |]E( e ^( a; i'- a; i)/^)E(e lr ' ia ' (y i'- %)/Z0 )| = 0, j V 3 

if ( |28| holds. In this case, 

< 2K 2 /^np 

with probability 1 — 26. 



□ 



5.2. Proof of Lemma [2l Lower bound. 

Proof. The Berry-Esseen theorem [25] states that the distribution of the sum of m inde- 
pendent and identically distributed zero-mean random variables normalized by its standard 

where a 2 and 



deviation, differs from the unit Gaussian distribution by at most Cp/{o" 
p are respectively the variance and the absolute third moment of the parent distribution, 
and C is a distribution-independent absolute constant which is not greater than 0.7655 [32J. 
We shall apply the Berry-Esseen theorem to the two summations, denoted by S p and T n 



respectively, on the right hand side of (52). 



The complex- valued random variables involved can be treated as M 2 -valued random vari- 



ables. Under (28) the variance of these random variables is 1/2 and the absolute third 
moment is 4/ (3n). 

Let Fi, F2 be the cumulative distributions of the real and imaginary parts of p~ 1 / 2 S p 
and Gi, G2 the cumulative distributions of the real and imaginary parts of rr x l 2 T n . Let 
be the cumulative distribution of the standard normal random variable. We have by the 
Berry-Esseen theorem 



(61) 
(62) 



sup !*;(«)-¥(*) I < f 8 ^, < = i,2 

sup \Gi(t)-9{t)\ < f 8 ^, < = 1,2. 



Since C < 0.7655, we can replace the right hand side of (61) and (62)) by p 1 ^ 2 and n 1 / 2 
respectively for the sake of notational simplicity. Hence 



\Fi(t) - Fi(-t)\ < \V{t)-V{-t)\ + 

\Gi(t)-Gi(-t)\ < \V(t)-y(-t)\ + 
Wt. For small t > we can bound the above expressions by 

\F(t)-F t (-t)\ < t ] /l + ^ 



\Gi{t)-Gi(-t)\ < t 
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n 



2 



n 



which imply 



and consequently 
(63) P 



P 
P 



S p T n \ 1t\ti 



P~ 1/2 \S P \ < ty/2 
n~ 1/2 \T n \ < tV2 



> \ 1-2U 



np yjnp 
which is what we want to prove. 



< 2t\ 

< 2t\ 

4 

71 yjp 



2 




4 




+ 




71 


Vp 




+ 


4 


7T 





1 - 2U 



n 



□ 



6. Average coherence bound: proof of Lemma [3] 



Proof. Write 



1 



max 

N-l j> 



n p 



, i = (k — l)n + I 



EEE^*< 

1=1 k=ljjLj> 

and consider the sums over k and j simultaneously with a fixed j' and fixed n sensor locations. 
This is a summation of p ■ N independent random variables 4>if4>ij each bounded by n~ 1 p~ 1 
in absolute value. Note that 

(64) E^[«/»*,0y]=O, 

since 6kj are uniformly distributed in [0, 2ir]. Applying Hoeffding inequality with 

c i 



(AT_ !)l/2p3/2 n ' 



c> 



we have 
(65) 



P. 



N — 1 



EE 

fc=l jVj' 



> 



Cl 



(N - l)l/2 p l/2 



n 



< 4e" Cl 



where Pg jr) is the probability conditioned on fixed £ = (£j),T) = (rjj) G M n . In analyzing the 
sum over 1 = 1 n we shall restrict to the event 



A 



= [9 kj ] : 



N-l 



< 



EE 

fc=i &j> 

Since there are at most N possible sensor locations, by (65) 
(66) 



npi/2(jv - 1)V2 

locatio) 
P(^l c ) < 4iW 



for almost all sensor locations. 



where A c denotes the complement of A. 
Let 



and E_4 is the expectation conditioned on the event A. 
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We proceed with the following estimate 



max 

r 



(Zj'l ~ ^A Z fl) 



1=1 



> 



C-2 



P 



A 



max 

i' 



^np(N - 1 



1=1 



> 



c 2 



A'- 



max 

r 



- ^A Z j'l, 



1=1 



y/np{N-l] 

C 2 

> 



F(A) 



(67) 



max 

3' 



( z j'i - ^A z j'i] 



i=i 



> 



^np(N - 1) 

C2 



+ 4Ne~ c \ ci,c 2 >0 



y/np(N-r 

by (66) where P^ and P^ c are respectively the probabilities conditioned on the events A and 
Applying Hoeffding's inequality with 

C 2 

t = pl/2(AT- l)l/2 n 3/2 

to estimate the first term on the right hand side of (67), we obtain 



P 



A 



(Zj'l - ^AZj'l 



1=1 



> 



c 2 



^Jnp(N - 1) 



< 4e - C |/(2 Cl ) : 



Maximizing over j' = 1, m and using the union bound we then arrive at 



(6£ 



P 



A 



max 
j' 



y (Zj'z - ^A z j'i] 



i=i 



> 



(-2 



Using (67) and (68) with 



C 2 = C4 



/JV- 1 
np 



c 2 = 2cj, c> 



we have 



P 



max 
j' 



> 



np 



_ c / N-l 

< 8Ne 2 V »p 



c> 0, 



which is what we set out to prove. 
Note that 

EqZm = — E ( e «6«-(ay-<'y)/*>) e (e 

TO V 7 V 

where Ee is the expectation conditioned on = (9kj) G C pxn . If 

1 Ai 



), Vj' = l,...,m, / = !,..., to, jVj 



p A2 



G N 



then 



E e 2j/ Z = 0, Vj' = 1, ...,rro, Z = 1,...,to 
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and hence 



^A z fi = 0, Vj' = l,...,m, 1 = 1, ...,n. 



□ 



7. Operator norm bound: proof of Lemma [4] 

Proof. It suffices to show that the matrix $ satisfies 

ii n P 



(69) 



N 



:$$*-I np || 2 <l 



where I np is the np x np identity matrix with the corresponding probability bound. Since 
the diagonal elements of ^2$$* are unity, (69) would in turn follow from 



(70) 



M (<&*)< 



1 



by the Gershgorin circle theorem. 
The pairwise coherence has the form 



N 



1 

N 



np — 1 



N 



3=1 



There are two cases: (i) k ^ k! , (ii) k = k',l I'. 

For case (i), 6^ — O^j are independent random variables for j = 1,...,N. Applying 
Hoeffding inequality to 



7V 



we obtain 
(71) 

Set t = a/y/N, we have 

P 



e *uxj (f ,/ )/z giw% -r)i)/z e %(Q ki -6 k ij) 

3=1 



P 



iV 



Zjvl > t 



< 4e 



-7Vt 2 



np 
~N 



N 



> 



a 



< 4e" a 



and thus 
(72) 



P 



sup 

Vij' 



TV 

E< 

3=1 



> 



a 



< 2n 2 p{p - l)e~ a 



by the union bound. 

For case (ii), 9^ — 6 k ij = and Z N becomes a geometric series 



Z N = ; ^ —, X 



I _ ^{^i-^t/zo 
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Thus, 



np 
iV 



N 






1 

< — 
- N 





sm 



2z 



sin 



2^o 



sm 



220 



sm 



2Zo 



Let 



k = mm mm 



Xz 



Km 1 - m) 



\z 



Clearly k is nonzero with probability one. For I ^ I' the probability density functions 
(PDF) for the random variables 



Az 



are either the symmetric triangular distribution or its self-convolution supported on [— 2p _1 , 2p _1 ] 
In either case, their PDFs are bounded by p. Hence the probability that {k > /?} for small 
(3 > is larger than 

(1 - 2pf3) n{n - 1)/2 > 1 - f3pn(n - 1) 

where the exponent counts the number of distinct unordered pairs (I, I'). Note that the above 
analysis is independent of k = k'. Since sin8 > #,V# £ [0, 7r/2] we have that 



< {3pn(n — 1). 







N 




(73) P 


np 
sup — 

h=k' 
_ l^l 1 


3=1 


TT 2 

~ 4N(3 2 



Setting 
(74) 



max 



a 



< 



71 



■s/N J 4Nf3 2 np-l 



and using (|72|) and (73|) we have 
np 



(75) 



P 



sup 



JV 
3=1 



> 



1 



np — 1 



< f3pn(n - 1) + 2n 2 pG» - l)e" a . 



As a consequence, 

AT 



p 



np 



hrfh 

3=1 



> 



np — 1 



< pn{n — 1)- 



7r / np — 1 

AT 



+ 2n 2 p{p — l)e ("p- 1 ) 



by maximizing the right hand side of ( 75 ) under the constraint ( 74 ) . 
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Figure 3. The Lasso performance comparison between RI with n = 11, p = 6 
and MR with n = 11. The vertical axis is for the success probability and 
the horizontal axis is for the number of objects. The success probability is 
estimated from 1000 independent trials. 



8. Numerical simulations 



We use two numerical settings: the diffraction-limited case when (28) is satisfied (Figure 
|5l [6]) and the under-resolved case when the ratio in (28) is smaller than unity (Figure 



7j. 

For the diffraction-limited case we set Zq = 10000 and A = 0.1 for the search domain 
[— 250, 250] 2 with t = 10. The targets are i.i.d. uniform random points in the grid with 
amplitudes in the range [1,2]. We randomly select sensor locations from [— 50, 50] 2 with the 
aperture A = 100 satisfying (28). With these parameters 



(A + £VN) 



1.3 



the condition (23) is barely satisfied. For the Lasso solution we have used the Matlab code 



Subspace Pursuit (available at http://igorcarron.googlepages.com/cscodes) 



We use the true Green function (20) in the computation of scattered waves and in recovery 



the exact Green function as well as its paraxial approximation to construct the sensing 
matrix (for comparison). In other words, we allow model mismatch between the forward 
and inversion steps. 
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Figure 4. The numbers of recoverable (by the Lasso) objects for RI with 
p = (n + l)/2 and MR as n varies. The curves indicate a quadratic behavior 
predicted by the theory. The difference between recoveries with the exact and 
paraxial Green functions is negligible in both the RI and MR set-ups. 



In the first set of simulations, we compare the performances of the Lasso for two imaging 
set-ups: one with random illumination (RI) and the other with multi-static responses (MR). 
As Figure [3] shows, the RI set-up has a higher success probability than the MR set-up. 
Another comparison is shown in Figure [4] in terms of the number of recoverable objects over 
a range of n. The quadratic behavior is consistent with the prediction of (39) and (36). The 
difference between the exact and paraxial Green functions recoveries is negligible in both the 
RI and MR set-ups. For a given n, the Lasso with the RI set-up recovers a higher number 
of objects than does the Lasso with the MR set-up. 

Figure [5] compares the performances of the Lasso (top panel) and OST (bottom panel) 
in terms of the number of recoverable objects for a fixed np = 600 but variable n. Clearly, 
the Lasso can recover far more objects exactly than does the OST. For a fixed np the 
performance for each method appears relatively constant over the whole range of n. For 
small n, the performance curves of both methods indicate superresolution. As noise level 
increases the Lasso performance decays (Figure [6]) . 

To further understand the superresolution effect of random illumination, we consider the 
set-up with zq = 25000, A = 0.4 for which the ratio in (28) is 0.1. This is an under-resolved 
case whose performance is shown in Figure [7j In contrast to the diffraction-limited case 
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Figure 5. The number of recoverable objects as a function of the number of 
sensors n = 1,2,3,4,5, 6, 8, 10, 12, 15, 20, 24, 25, 30, 40, 50, 60, 75, 100, 120, 
150, 200, 300, 600 with np = 600 fixed. The top panel is for the Lasso and the 
bottom panel for OST. The left ends of both curves indicate superresolution. 



(Figure [5]), the number of recoverable objects in the under-resolved case decays rapidly as 
p decreases (n increases). To maintain high performance in the under- resolved case, it is 
necessary that p^> 1. The number of recoverable objects is calculated based on 90% recovery 
of 100 independent trials. 

We demonstrate in Figures [8- VL the performance for extended objects in the presence of 



^=K + ^ 2 )^, p = 5%,20% 
V2 y/np 



external noise of the form 



where p is the percentage of noise in each entry of the data vector and V\, z/ 2 are i.i.d. uniform 
random variables in [0, 1]. 

Figure M shows the original 40 x 80 pixel image (left) and its reconstructions (middle panel, 



5% noise; right panel, 20% noise) by the BPDN solver YALL1 (http : //yalll . blogs . rice . edu/) 



using one sensor and 500 random illuminations while Figure [9] shows the results with one 
illuminations and 500 randomly distributed sensors. 



Figure 10 shows the original 70 x 70 pixel image (left), the Shepp-Logan phantom, and its 



reconstructions (middle panel, 5% noise; right panel, 20% noise) by the total-variation min- 



imization [TH [3TJ solver TVAL3 (http://www.caam.rice.edu/~optimization/Ll/TVAL3/) 



using one sensor and 1000 random illuminations while Figure [TT] shows the results with one 
illumination and 1000 randomly distributed sensors. 
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Figure 6. Noisy recovery by the Lasso for n = 1, 30, 100, 600 and with np = 
600 fixed. The noise is given by the circularly random Gaussian noise of 
magnitude cr||y||2 where a is the horizontal coordinate. Note that in this case 
E\\E\\ 



npa 2 \\Y\\l. 



The low pixel numbers are chosen to reduce the run time of the programs. 



For the one-illumination reconstructions (Figures P and 11), the classical resolution crite- 



rion (28) is met. Note, however, that the Shepp-Logan phantom is not in the class of sparse 
extended objects analyzed in Section [|] because the object support covers more than 50% of 
the domain (only the gradient is sparse). As a result, the same percentage of noise represents 
a greater amount of noise in the case of Shepp-Logan phantom and has a more serious effect 



on performance (Figures 10 and Oil right panels). 



9. Conclusion 

We have proposed a new approach to superresolving point and extended objects based on 
random illumination and compressed sensing reconstruction. 

We have proved that in the diffraction-limited case both the Lasso and the OST with 
random illumination can exactly localize s = 0(m) objects where the number of data m is 
the product of the numbers of random probes and sensors. For the under-resolved case where 
the Rayleigh resolution limit is broken, the Lasso still has a similar performance guarantee 
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Figure 7. The number of recoverable objects in the under-resolved 

a function of the number of sensors n = 1,2,3,4,5, 6, 8, 10, 12, 15, 20, 24, 25, 

30, 40, 50, 60, 75, 100, 120, 150, 200, 300, 600 with np = 600 fixed. 





LiC 



Figure 8. The original 40 x 80 pixel image (left) and the BPDN reconstruc- 
tions (middle panel, 5% noise; right panel 20% noise) with one sensor and 500 
random illuminations. 



if the number of random illuminations is sufficiently large. It is possible to extend the OST 
result to the under-resolved case which is omitted here to simplify the presentation. 

Numerical evidence supports our theoretical prediction and confirms the superiority of the 
Lasso to the OST in the set-up with random illumination. 

We have also shown that the BPDN is suitable for imaging extended objects and have 
provided numerical examples to demonstrate its performance. 
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Figure 9. The original 40 x 80 pixel image (left) and the BPDN reconstruc- 
tions (middle panel, 5% noise; right panel, 20% noise) with one illumination 
and 500 randomly distributed sensors. 
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Figure 10. The original 70 x 70 pixel image (left), the Shepp-Logan phantom, 
and the TV-minimization reconstructions (middle panel, 5% noise; right panel, 
20% noise) with one sensor and 1000 random illuminations. 
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Figure 1 1 . The original 70 x 70 pixel image (left) , the Shepp-Logan phantom, 
and the TV-minimization reconstructions (middle panel, 5% noise; right panel, 
20% noise) with one illumination and 1000 randomly distributed sensors. 




The superresolution effect with random illumination revealed here contrasts with the sub- 
wavelength resolution with deterministic near-field illumination studied in [20J. 
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Finally we note that in our approach it is essential to measure the wave field. For intensity- 
only measurements, additional techniques such as interferometry or phase retrieval methods 
are necessary for object reconstruction. 



Acknowledgement I am grateful to Mike Yan for producing Figures [3](7] and Hsiao- 
Chieh Tseng for producing Figures |8]fTl~1 of Section [8j 
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